Code
library(readxl)
library(EpiEstim)
library(ggplot2)
library(dplyr)
library(janitor)
library(tidyr)
library(knitr)
library(lubridate)
df1 <- read_xlsx("data/data.xlsx", sheet = "data")
df1 <- as.data.frame(df1)library(readxl)
library(EpiEstim)
library(ggplot2)
library(dplyr)
library(janitor)
library(tidyr)
library(knitr)
library(lubridate)
df1 <- read_xlsx("data/data.xlsx", sheet = "data")
df1 <- as.data.frame(df1)Clean data
df1 <- df1 %>% clean_names()
df <- df1 %>% rename(dates = ngay_khoi_phat,
ngaysinh = ngay_sinh,
ngaynv = ngay_nhap_vien_kham,
tiemchung = tinh_trang_tiem_chung,
gioi = gioi_tinh)
sum(is.na(df$dates))[1] 173
# Lấy ngày nhập viện thay cho ngày khởi phát của các ca missing
df$dates[is.na(df$dates)] <- df$ngaynv[is.na(df$dates)]
sum(is.na(df$dates))[1] 0
df$cd <- ifelse(df$phan_loai_chan_doan == "Loại trừ sởi", NA, df$phan_loai_chan_doan)
df <- df[,c("dates", "ngaysinh", "tiemchung", "cd")]
# Loại những ca Sởi loại trừ
df <- na.omit(df)
df_convert <- df %>% group_by(dates) %>%
summarise(I = n())
df_complete <- df_convert %>%
complete(dates = seq(min(dates), max(dates), by = "day")) %>%
replace_na(list(I = 0))
df_complete$dates <- as.Date(df_complete$dates)Ước tính hệ số lây nhiễm Rt
mod_raw <- estimate_R(
incid = df_complete, # Dữ liệu từ 04/03/2024 - end
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25
)
)
)Default config will estimate R on weekly sliding windows.
To change this change the t_start and t_end arguments.
Warning in estimate_R_func(incid = incid, method = method, si_sample = si_sample, : You're estimating R too early in the epidemic to get the desired
posterior CV.
df_rt_raw <- mod_raw$R
df_rt_raw$dates <- mod_raw$dates[df_rt_raw$t_end]
df_rt_raw$q1_rt <- df_rt_raw$`Quantile.0.025(R)`
df_rt_raw$q3_rt <- df_rt_raw$`Quantile.0.975(R)`
df_rt_raw$rt <- df_rt_raw$`Mean(R)`library(patchwork)
p_hist_raw <- ggplot(df_complete, aes(x = dates, y = I)) +
geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
labs(x = "Day", y = "Incidence") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-03-04"), ymd("2024-09-14")))
p_raw <- ggplot(df_rt_raw, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "04/03/2024 - 14/09/2024") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-03-04"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 20))
p_hist_raw/p_raw# Lọc dữ liệu từ 23/05/2024 - end
df_filter <- filter(df_complete, dates >= "2024-05-30")mod_1w <- estimate_R(
incid = df_filter,
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25
)
)
)Default config will estimate R on weekly sliding windows.
To change this change the t_start and t_end arguments.
# Lấy dữ liệu ra để vẽ Rt
df_rt_1w <- mod_1w$R
df_rt_1w$dates <- mod_1w$dates[df_rt_1w$t_end]
df_rt_1w$q1_rt <- df_rt_1w$`Quantile.0.025(R)`
df_rt_1w$q3_rt <- df_rt_1w$`Quantile.0.975(R)`
df_rt_1w$rt <- df_rt_1w$`Mean(R)`library(patchwork)
# Biểu đồ đường cong dịch
p_hist <- ggplot(df_filter, aes(x = dates, y = I)) +
geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
labs(x = "Day", y = "Incidence") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14")))
# Biểu đồ cho 1 tuần
p_1w <- ggplot(df_rt_1w, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "30/05/2024 - 14/09/2024") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 20))library(patchwork)
(p_hist_raw + p_hist) / (p_raw + p_1w)esimate_R, chúng ta có thể tùy chọn Sliding window (SW) khác nhau để thể hiện kết quả ước tính Rt. Sliding window là khoảng thời gian chúng ta dùng để xác định Rt, ví dụ nếu chúng ta chọn SW = 7 thì giá trị Rt sẽ được tính bằng trung bình của 7 ngày. Vậy chúng ta nên chọn *SW bằng bao nhiêu? Chúng ta vẽ thử 4 trường hợp của SW lần lượt là: mỗi ngày, mỗi tuần, mỗi 2 tuần và mỗi tháng để lựa chọn.# SW = 1
t_start <- seq(2, nrow(df_filter)-1)
t_end <- t_start + 1
mod_day <- estimate_R(
incid = df_filter,
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25,
t_start = t_start,
t_end = t_end
)
)
)
# Lấy dữ liệu ra để vẽ Rt
df_rt_day <- mod_day$R
df_rt_day$dates <- mod_day$dates[df_rt_day$t_end]
df_rt_day$q1_rt <- df_rt_day$`Quantile.0.025(R)`
df_rt_day$q3_rt <- df_rt_day$`Quantile.0.975(R)`
df_rt_day$rt <- df_rt_day$`Mean(R)`# SW = 14
t_start <- seq(2, nrow(df_filter)-13)
t_end <- t_start + 13
mod_2w <- estimate_R(
incid = df_filter,
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25,
t_start = t_start,
t_end = t_end
)
)
)
# Lấy dữ liệu ra để vẽ Rt
df_rt_2w <- mod_2w$R
df_rt_2w$dates <- mod_2w$dates[df_rt_2w$t_end]
df_rt_2w$q1_rt <- df_rt_2w$`Quantile.0.025(R)`
df_rt_2w$q3_rt <- df_rt_2w$`Quantile.0.975(R)`
df_rt_2w$rt <- df_rt_2w$`Mean(R)`# SW = 28
t_start <- seq(2, nrow(df_filter)-27)
t_end <- t_start + 27
mod_month <- estimate_R(
incid = df_filter,
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25,
t_start = t_start,
t_end = t_end
)
)
)
# Lấy dữ liệu ra để vẽ Rt
df_rt_month <- mod_month$R
df_rt_month$dates <- mod_month$dates[df_rt_month$t_end]
df_rt_month$q1_rt <- df_rt_month$`Quantile.0.025(R)`
df_rt_month$q3_rt <- df_rt_month$`Quantile.0.975(R)`
df_rt_month$rt <- df_rt_month$`Mean(R)`library(patchwork)
# Biểu đồ đường cong dịch
p_hist <- ggplot(df_filter, aes(x = dates, y = I)) +
geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
labs(x = "Day", y = "Incidence") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14")))
# Biểu đồ cho 1 tuần
p_1w <- ggplot(df_rt_1w, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "1 Week") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 8))
# Biểu đồ cho 1 ngày
p_day <- ggplot(df_rt_day, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "1 Day") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 10))
# Biểu đồ cho 2 tuần
p_2w <- ggplot(df_rt_2w, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "2 Week") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 10))
# Biểu đồ cho 1 tháng
p_month <- ggplot(df_rt_month, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "1 Month") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 5))
(p_hist + p_hist)/(p_day + p_1w)library(patchwork)
(p_hist + p_hist)/(p_2w + p_month)Rt cao khi số ca bắt đầu tăng đột ngột và liên tục, càng về sau Rt giảm xuống gần bằng 1 nhưng không duy trì được giá trị Rt < 1 chứng tỏ dịch vẫn đang diễn ra. Với SW = 1 tuần thì Rt ~ 14; Sliding window = 2 tuần thì Rt ~ 19; các khoảng Sliding window khác có giá trị Rt thấp hơn rất nhiều vào ngày ước tính đầu tiên. Theo lý thuyết, Rt sẽ cao khoảng thời gian đầu do dịch mới bùng phát mà 2 tuần trước đó không có ca bệnh xuất hiện liên tục. Do đó, chúng ta chọn giữa SW = 7 hoặc SW = 14 để ước tính.
Tuy nhiên, kết quả cho thấy Rt được ước tính với SW là 2 tuần cho kết quả ổn định hơn, giá trị không lên xuống thất thường nên phù hợp để xem xét và đánh giá diễn tích của dịch. Ngoài ra, khoảng tin cậy của Rt trong thời đầu là quá lớn (khi SW là 1 tuần) nên sẽ không đảm bảo tính tin cậy của Rt.
Ước tính Rt với sliding window là 2 tuần
library(patchwork)
p_hist <- ggplot(df_filter, aes(x = dates, y = I)) +
geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
labs(x = "Day", y = "Incidence") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14")))
p_2w <- ggplot(df_rt_2w, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "13/06 - 14/09") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(breaks = seq(0,10,by = 2))
p_hist/p_2wlibrary(knitr)
kable(df_rt_2w[,c("dates", "rt")])| dates | rt |
|---|---|
| 2024-06-13 | 7.0698366 |
| 2024-06-14 | 5.6690111 |
| 2024-06-15 | 4.3508564 |
| 2024-06-16 | 3.2821600 |
| 2024-06-17 | 2.4718480 |
| 2024-06-18 | 2.1445391 |
| 2024-06-19 | 1.8320200 |
| 2024-06-20 | 1.7193623 |
| 2024-06-21 | 1.6256121 |
| 2024-06-22 | 1.5422473 |
| 2024-06-23 | 1.5176277 |
| 2024-06-24 | 1.4954698 |
| 2024-06-25 | 1.3944655 |
| 2024-06-26 | 1.4405630 |
| 2024-06-27 | 1.3418515 |
| 2024-06-28 | 1.1143724 |
| 2024-06-29 | 1.0995479 |
| 2024-06-30 | 1.2085188 |
| 2024-07-01 | 1.4736674 |
| 2024-07-02 | 1.5490903 |
| 2024-07-03 | 1.5791076 |
| 2024-07-04 | 1.5973473 |
| 2024-07-05 | 1.5671345 |
| 2024-07-06 | 1.4948245 |
| 2024-07-07 | 1.4237618 |
| 2024-07-08 | 1.3243248 |
| 2024-07-09 | 1.1656958 |
| 2024-07-10 | 1.0137158 |
| 2024-07-11 | 1.0568080 |
| 2024-07-12 | 1.0603317 |
| 2024-07-13 | 1.1164884 |
| 2024-07-14 | 1.0754648 |
| 2024-07-15 | 0.8953696 |
| 2024-07-16 | 0.9462627 |
| 2024-07-17 | 1.0530632 |
| 2024-07-18 | 1.0882127 |
| 2024-07-19 | 1.1613209 |
| 2024-07-20 | 1.2219629 |
| 2024-07-21 | 1.2947699 |
| 2024-07-22 | 1.4031893 |
| 2024-07-23 | 1.4613128 |
| 2024-07-24 | 1.5730353 |
| 2024-07-25 | 1.6793297 |
| 2024-07-26 | 1.7178643 |
| 2024-07-27 | 1.6554286 |
| 2024-07-28 | 1.6616364 |
| 2024-07-29 | 1.7875178 |
| 2024-07-30 | 1.7556702 |
| 2024-07-31 | 1.6092751 |
| 2024-08-01 | 1.4647907 |
| 2024-08-02 | 1.3967257 |
| 2024-08-03 | 1.3536969 |
| 2024-08-04 | 1.3559618 |
| 2024-08-05 | 1.3996821 |
| 2024-08-06 | 1.5009996 |
| 2024-08-07 | 1.6549297 |
| 2024-08-08 | 1.7120074 |
| 2024-08-09 | 1.8062972 |
| 2024-08-10 | 1.8852197 |
| 2024-08-11 | 1.9690807 |
| 2024-08-12 | 1.8912403 |
| 2024-08-13 | 1.9860504 |
| 2024-08-14 | 1.9307622 |
| 2024-08-15 | 1.9485278 |
| 2024-08-16 | 1.9479556 |
| 2024-08-17 | 1.9839644 |
| 2024-08-18 | 1.9756048 |
| 2024-08-19 | 1.9431634 |
| 2024-08-20 | 1.7682855 |
| 2024-08-21 | 1.6166982 |
| 2024-08-22 | 1.5031903 |
| 2024-08-23 | 1.4052603 |
| 2024-08-24 | 1.3527849 |
| 2024-08-25 | 1.3302986 |
| 2024-08-26 | 1.3423795 |
| 2024-08-27 | 1.2956510 |
| 2024-08-28 | 1.3565252 |
| 2024-08-29 | 1.3408987 |
| 2024-08-30 | 1.2963237 |
| 2024-08-31 | 1.2655504 |
| 2024-09-01 | 1.2486194 |
| 2024-09-02 | 1.2373010 |
| 2024-09-03 | 1.3128349 |
| 2024-09-04 | 1.3294864 |
| 2024-09-05 | 1.3233001 |
| 2024-09-06 | 1.2926711 |
| 2024-09-07 | 1.2373607 |
| 2024-09-08 | 1.1799842 |
| 2024-09-09 | 1.1565229 |
| 2024-09-10 | 1.1126126 |
| 2024-09-11 | 1.0638421 |
| 2024-09-12 | 1.0497415 |
| 2024-09-13 | 1.0363817 |
| 2024-09-14 | 1.0485225 |